#include "mcmcParticleSurface.h"


/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
// construction & destruction
MCMC_ControlPoints::MCMC_ControlPoints()
{
	pmx_A = NULL;
}

MCMC_ControlPoints::~MCMC_ControlPoints()
{
	if (pmx_A != NULL)
		delete pmx_A;
}

MCMC_ControlPoints& MCMC_ControlPoints::operator=(const MCMC_ControlPoints& a)
{
	if (this != &a)
	{
		nrow			= a.nrow;
		ncol			= a.ncol;
		mx_P			= a.mx_P;
		vc_u			= a.vc_u;
		vc_w			= a.vc_w;
		mx_par			= a.mx_par;
		num_pixels		= a.num_pixels;
		belongated		= a.belongated;
		pmx_medial_row	= a.pmx_medial_row;
		pmx_medial_col	= a.pmx_medial_col;
		pmx_medial_dist = a.pmx_medial_dist;
		vc_a			= a.vc_a;
		vc_c			= a.vc_c;
	}

	return *this;
}

int MCMC_ControlPoints::TotalCodingBits(void)
{
	int bits;

//	bits = 2*sizeof(byte)*vc_u.size()*vc_w.size();	// u and w matrices
	
	bits = sizeof(float)*mx_par.size();	// u and w matrices

	return bits;
}

// operations
void MCMC_ControlPoints::Initialize(mcmcMatrix<int> &mx_label, const int label, mcmcMatrix<int> *pmx_taken, Raster<float> &image,
									irmm_Rect *pbound)
{
	int			i,j;
	irmm_Rect	rect_bound;

	nrow = mx_label.rows();
	ncol = mx_label.cols();

	// how many number of pixels belong to this cluster
	num_pixels = 0;

	if (pbound != NULL)
		rect_bound = *pbound;
	else
		rect_bound.Set(0, 0, mx_label.cols()-1, mx_label.rows()-1);

	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_label(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
				num_pixels++;
	if (num_pixels < 10 )
		return;

	belongated		= false;
	pmx_medial_row	= NULL;
	pmx_medial_col	= NULL;
	pmx_medial_dist = NULL;

	mcmcMatrix<double>	A,AT,P,B;
	double			x,y;
	int			t;

	A.SetDimension(num_pixels, 6);
	AT.SetDimension(6, num_pixels);
	P.SetDimension(num_pixels, 1);
	t = 0;
	
	// get each pixel one by one
	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_label(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
			{
				x = (double)j;  y = (double)i;
				A(t, 0) = x*x;  AT(0, t) = A(t, 0);
				A(t, 1) = x*y;	AT(1, t) = A(t, 1);
				A(t, 2) = y*y;  AT(2, t) = A(t, 2);
				A(t, 3) = x;    AT(3, t) = A(t, 3);
				A(t, 4) = y;    AT(4, t) = A(t, 4);
				A(t, 5) = 1.0;  AT(5, t) = A(t, 5);
							
				P(t) = image(j, i);
				t++;
			}

	B = AT*A;
	
	mx_par = B.Inv()*AT*P;
}

void MCMC_ControlPoints::Initialize(LinkList<PMCMC_PixelSet> &list_domain, mcmcMatrix<int> *pmx_taken, Raster<float> &image)
{
	int			i,j;
	irmm_Rect	rect_bound;
	bool		bgo;

	nrow = image.rows();
	ncol = image.cols();

	// how many number of pixels belong to this cluster
	num_pixels = 0;

	bool			bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	if (pmx_taken == NULL)
	{
		bgo_domain = list_domain.GetFirst(ppixel_set);
		while (bgo_domain)
		{
			if (ppixel_set != NULL)
				num_pixels += ppixel_set->Num();	
			
			bgo_domain = list_domain.GetNext(ppixel_set);
		}
		
	}
	else
	{
		bgo_domain = list_domain.GetFirst(ppixel_set);
		while (bgo_domain)
		{
			if (ppixel_set != NULL)
			{			
				bgo = ppixel_set->GetFirst(j, i);
				while (bgo)
				{
					if (pmx_taken->Valid(j, i) && pmx_taken->Data(j,i)==0)
						num_pixels++;
					bgo = ppixel_set->GetNext(j, i);
				}
			}
			bgo_domain = list_domain.GetNext(ppixel_set);
		}
	}

	belongated		= false;
	pmx_medial_row	= NULL;
	pmx_medial_col	= NULL;
	pmx_medial_dist = NULL;

	if (num_pixels < 10 )
		return;
	
	mcmcMatrix<double>	A,AT,P,B;
	double			x,y;
	int			t;

	A.SetDimension(num_pixels, 6);
	AT.SetDimension(6, num_pixels);
	P.SetDimension(num_pixels, 1);
	t = 0;
	
	bgo_domain = list_domain.GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{			
			bgo = ppixel_set->GetFirst(j, i);
			// get each pixel one by one
			while (bgo)
			{
				if (pmx_taken==NULL || (pmx_taken->Valid(j, i) && pmx_taken->Data(j,i)==0))
				{
					x = (double)j;  y = (double)i;
					A(t, 0) = x*x;  AT(0, t) = A(t, 0);
					A(t, 1) = x*y;	AT(1, t) = A(t, 1);
					A(t, 2) = y*y;  AT(2, t) = A(t, 2);
					A(t, 3) = x;    AT(3, t) = A(t, 3);
					A(t, 4) = y;    AT(4, t) = A(t, 4);
					A(t, 5) = 1.0;  AT(5, t) = A(t, 5);
									
					P(t) = image(j, i);
					t++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = list_domain.GetNext(ppixel_set);
	}

	B = AT*A;
	
	mx_par = B.Inv()*AT*P;
}


void MCMC_ControlPoints::Initialize(mcmcMatrix<int> &mx_cluster1, const int label, mcmcMatrix<int> *pmx_taken,
								   Raster<float> &raster1, mcmcMatrix<int> &mx_medial_row, mcmcMatrix<int> &mx_medial_col,
								   mcmcMatrix<int> &mx_medial_dist, irmm_Rect *pbound)
{
	int			i,j;
	irmm_Rect	rect_bound;

	nrow = mx_cluster1.rows();
	ncol = mx_cluster1.cols();

	pmx_medial_row = &mx_medial_row;
	pmx_medial_col = &mx_medial_col;
	pmx_medial_dist = &mx_medial_dist;

	// how many number of pixels belong to this cluster
	num_pixels = 0;

	if (pbound != NULL)
		rect_bound = *pbound;
	else
		rect_bound.Set(0, 0, mx_cluster1.cols()-1, mx_cluster1.rows()-1);

	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_cluster1(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
			{
				if (mx_medial_row(j,i)==j && mx_medial_col(j,i)==i)
					num_pixels++;
			}
	if (num_pixels < 7)
		return;

	belongated = true;

	mcmcMatrix<double>	A,AT,P,B;
	double			x,y,d;
	int			t;
	int				row,col,k,nb_row[4]={-1,0,0,1},nb_col[4]={0,-1,1,0};

	A.SetDimension(num_pixels, 6);
	AT.SetDimension(6, num_pixels);
	P.SetDimension(num_pixels, 1);
	t = 0;
	
	//roughly fit the surface
	// get each pixel one by one
	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_cluster1(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
			{
				if (mx_medial_row(j,i)==j && mx_medial_col(j,i)==i)
				{
					x = j;	y = i;
					
					A(t, 0) = x*x;  AT(0, t) = A(t, 0);
					A(t, 1) = x*y;	AT(1, t) = A(t, 1);
					A(t, 2) = y*y;  AT(2, t) = A(t, 2);
					A(t, 3) = x;    AT(3, t) = A(t, 3);
					A(t, 4) = y;    AT(4, t) = A(t, 4);
					A(t, 5) = 1.0;  AT(5, t) = A(t, 5);
					
					P(t) = raster1(j,i);
					t++;
				}
			}
	B = AT*A;
	
	mx_par = B.Inv()*AT*P;
	
	num_pixels = 0;
	// robust fitting
	double dvalue;
	int		num_medial=0;
	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_cluster1(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
			{
				x = j;	y = i;
				dvalue = mx_par(0)*x*x + mx_par(1)*x*y + mx_par(2)*y*y +
			  		     mx_par(3)*x   + mx_par(4)*y   +
				         mx_par(5);
/*				if (fabs(raster1(j,i)-dvalue)<10.0)
				{
					mx_medial_row(j, i) = j;
					mx_medial_col(j, i) = i;
				}
*/				if (mx_medial_row(j,i)==j && mx_medial_col(j,i)==i)
					num_pixels++;
			}
	mcmcVector<double> vc_num_temp;

	vc_a.SetRowVector(num_pixels);
	vc_a.InitValue(0.0);
	vc_c.SetRowVector(num_pixels);

	vc_num_temp.SetRowVector(vc_a.size());
	vc_num_temp.InitValue(0.0);

	k = 0;
	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_cluster1(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
			{
				if (mx_medial_row(j,i)==j && mx_medial_col(j,i)==i)
				{
					mx_medial_dist(j,i) = k;
					vc_c(k) = raster1(j,i);

					k++;
				}
			}

	for (j=rect_bound.top; j<=rect_bound.bottom; j++)
		for (i=rect_bound.left; i<=rect_bound.right; i++)
			if (mx_cluster1(j,i)==label && (pmx_taken==NULL || pmx_taken->Data(j,i)==0))
			{
				row = max(0, mx_medial_row(j,i));
				col = max(0, mx_medial_col(j,i));
				if (row!=j || col!=i)
				{
					k = mx_medial_dist(row, col);
					if (vc_a.Valid(k))
					{
						d = sqrt((float)((row-j)*(row-j) + (col-i)*(col-i)));
						vc_a(k) += raster1(j, i)-raster1(row, col);

						vc_num_temp(k) += d;
					}
				}
			}
	for (k=0; k<vc_num_temp.size(); k++)
	{
		if (vc_num_temp(k) > 0.0)
			vc_a(k) = vc_a(k)/vc_num_temp(k);
	}
}


void MCMC_ControlPoints::ComputeUW(const int j, const int i)
{
	if (num_pixels <= 0)
		return;				// no valid pixels exist

	double	u,w;

	u = ((double)j)/(nrow-1);
	w = ((double)i)/(ncol-1);
	
	vc_u(0) = pow(1-u,3.0);			vc_u(1) = 3.0*u*pow(1-u,2.0);
	vc_u(2) = 3.0*pow(u,2.0)*(1-u); vc_u(3) = pow(u,3.0);

	vc_w(0) = pow(1-w,3.0);			vc_w(1) = 3.0*w*pow(1-w,2.0);
	vc_w(2) = 3.0*pow(w,2.0)*(1-w); vc_w(3) = pow(w,3.0);
}

double	MCMC_ControlPoints::SurfaceValue(const int row, const int col)
{
	if (num_pixels <= 0)
		return 0.0;			// no valid pixels exist

/*
	double	value=0.0;
	int		m,n;

	ComputeUW(row,col);
	for (m=0; m<vc_u.size(); m++)
		for (n=0; n<vc_w.size(); n++)
			value += vc_u(m)*vc_w(n)*mx_P(m*vc_u.size()+n);
*/
	double x,y;
	
	if (!belongated)
	{
		x = (double)row;
		y = (double)col;

		if (mx_par.size() > 0)
			return mx_par(0)*x*x + mx_par(1)*x*y + mx_par(2)*y*y +
				   mx_par(3)*x   + mx_par(4)*y   +
				   mx_par(5);
		else
			return 0;
	}
	else
	{
		int		i,j,k;
		double	d,dvalue;

		j = max(0, pmx_medial_row->Data(row, col));
		i = max(0, pmx_medial_col->Data(row, col));
		
		k = pmx_medial_dist->Data(j, i);
		x = (double)row;
		y = (double)col;
		d = sqrt((float)((j-row)*(j-row)+(i-col)*(i-col)));
		
		if (mx_par.size() > 0)
		{
			dvalue = mx_par(0)*x*x + mx_par(1)*x*y + mx_par(2)*y*y +
					 mx_par(3)*x   + mx_par(4)*y   +
				     mx_par(5);
			if (vc_a.Valid(k))
				dvalue = vc_c(k)+vc_a(k)*d; 
			return dvalue;
		}
		else
			return 0;
	}
}



void MCMC_ControlPoints::Estimate(void)
{
	if (num_pixels <= 0)
		return;				// no valid pixels exist

/*
	mcmcMatrix<double> mx_temp,mx_AT;
	
	mx_AT = pmx_A->T();
	mx_temp = mx_AT*(*pmx_A);
	mx_P = mx_temp.Inv()*mx_AT*mx_B;

	if (pmx_A != NULL)
		delete pmx_A;
	pmx_A = NULL;
*/
}

/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
MCMC_Particle_Surface::MCMC_Particle_Surface()
{
	p_surface_model = NULL;
	p_synthesized   = NULL;
}

MCMC_Particle_Surface::~MCMC_Particle_Surface()
{
	if (p_surface_model != NULL)
		delete []p_surface_model;
	p_surface_model = NULL;

	if (p_synthesized != NULL)
		delete []p_synthesized;
	p_synthesized = NULL;
}

void MCMC_Particle_Surface::FitSurface(void)
{
	if (p_surface_model != NULL)
		delete []p_surface_model;
	p_surface_model = NULL;
	p_surface_model = new MCMC_ControlPoints[num_particle];
	// free synthesized images
	if (p_synthesized != NULL)
		delete []p_synthesized;
	p_synthesized = NULL;
	// allocate space
	p_synthesized = new Raster<float>[num_particle];

	for (int k=0; k<num_particle; k++)
	{
		p_surface_model[k].Initialize(mx_map_temp,k,NULL,prespFilters->ras_intensity);
		p_surface_model[k].Estimate();
		
		p_synthesized[k].SetDimension(prespFilters->rows(),prespFilters->cols());
		for (int j=0; j<prespFilters->rows(); j++)
			for (int i=0; i<prespFilters->cols(); i++)
				p_synthesized[k](j,i) = (float)p_surface_model[k].SurfaceValue(j,i);
	}
}


void MCMC_Particle_Surface::Set(MCMC_FilterResp &respFilters, mcmcMatrix<int> &clst, char *pfile_name)
{
	// initialization
	prespFilters = &respFilters;

	mx_map_temp = clst;
	num_particle = mx_map_temp.GetMax()+1;
	
	if (p_surface_model != NULL)
		delete []p_surface_model;
	p_surface_model = new MCMC_ControlPoints[num_particle];
	
	FitSurface();

	EM();

	CreateSalincyMap();
}

void MCMC_Particle_Surface::EM(void)
{
	int		i,j,k;
	double	dtemp;
	RF		rf;

	rf.Set(0,num_particle-1,1);

	for (int step=0; step<3; step++)
	{
		for (j=0; j<prespFilters->rows(); j++)
			for (i=0; i<prespFilters->cols(); i++)
			{
				for (k=0; k<num_particle; k++)
				{
					dtemp = DistanceToAParticle(prespFilters->ImageIntensity(j,i),k,j,i);
					if (dtemp>1e-6)
						rf.SetOneBin(k,1/dtemp);
					else
						rf.SetOneBin(k,1e6);
				}
				k = (int)rf.SamplingNonSetCdf();
				mx_map_temp(j,i) = k;
			}
		// setp 2:
		FitSurface();
	}

	// any two clusters are overlapping?
	mcmcVector<int> vc_sign(num_particle);
	vc_sign.InitValue(1);
	for (i=0; i<num_particle; i++)
	{
		if (vc_sign(i)==1)
		{
			for (j=i+1; j<num_particle; j++)
			{
				if (vc_sign(j)==1)
				{
					if (DistanceBtwnParticles(i,j)<10.0)
						vc_sign(j) = 0;		// tow clusters are too close
				}
			}
		}
	}
	mcmcMatrix<int> vc_index(num_particle,1);
	double		dist_temp;
	int			ipar_temp,m;

	m = 0;
	for (i=0; i<num_particle; i++)
		if (vc_sign(i)==1)
			vc_index(i) = m++;
	
	// relabel
	for (j=0; j<prespFilters->rows(); j++)
		for (i=0; i<prespFilters->cols(); i++)
		{
			dist_temp = 1e90;
			ipar_temp = -1;

			for (k=0; k<num_particle; k++)
			{
				if (vc_sign(k)==1)
				{
					dtemp = DistanceToAParticle(prespFilters->ImageIntensity(j,i),k,j,i);
					if (dtemp<dist_temp)
					{
						ipar_temp = k;
						dist_temp = dtemp;
					}
				}
			}
			assert(ipar_temp>=0);
			mx_map_temp(j,i) = vc_index(ipar_temp);
		}

	num_particle = m;
	FitSurface();

	Raster<uchar> aa;
	char		  name[200];
	aa.SetDimension(prespFilters->rows(), prespFilters->cols());
	for (k=0; k<num_particle; k++)
	{
		for (j=0; j<prespFilters->rows(); j++)
			for (i=0; i<prespFilters->cols(); i++)
				aa(j,i) = (uchar)(p_synthesized[k](j,i));
		sprintf(name,"a_%d.ppm",k);
		WriteToPPM(aa,name);
	}
}

double MCMC_Particle_Surface::DistanceToAParticle(double intensity, const int ipar, const int j, const int i)
{
	return fabs(p_synthesized[ipar](j,i)-intensity);
}

double MCMC_Particle_Surface::DistanceBtwnParticles(const int ipar1, const int ipar2)
{
	double	dis1,dis2;
	int	num1,num2;

	//initialization
	dis1 = 0.0;
	dis2 = 0.0;
	num1 = 0;
	num2 = 0;

	for (int j=0; j<prespFilters->rows(); j++)
		for (int i=0; i<prespFilters->cols(); i++)
		{
			if (mx_map_temp(j,i)==ipar1)
			{
				dis1 += fabs(p_synthesized[ipar1](j,i)-p_synthesized[ipar2](j,i));
				num1++;
			}
			else if (mx_map_temp(j,i)==ipar2)
			{
				dis2 += fabs(p_synthesized[ipar1](j,i)-p_synthesized[ipar2](j,i));
				num2++;
			}
		}
	if (num1 <= 0)
		dis1 = 1e5;
	else
		dis1 = (dis1/(double)num1);
	if (num2 <= 0)
		dis2 = 1e5;
	else
		dis2 = (dis2/(double)num2);
	return min(dis1,dis2);
}


void MCMC_Particle_Surface::CreateSalincyMap(void)
{
	RF		rf;
	int		i,j,k;
	double	dtemp;

	if (pmx_saliency_map != NULL)
		delete []pmx_saliency_map;
	pmx_saliency_map =  new mcmcMatrix<float>[num_particle];
	for (k=0; k<num_particle; k++)
		pmx_saliency_map[k].SetDimension(prespFilters->rows(),prespFilters->cols());

	rf.Set(0,num_particle-1,1);
	
	for (j=0; j<prespFilters->rows(); j++)
		for (i=0; i<prespFilters->cols(); i++)
			{
				for (k=0; k<num_particle; k++)
				{
					dtemp = DistanceToAParticle(prespFilters->ImageIntensity(j,i),k,j,i);
					if (dtemp>1e-6)
						rf.SetOneBin(k,1/dtemp);
					else
						rf.SetOneBin(k,1e6);
				}
				rf.normalize();
				for (k=0; k<num_particle; k++)
					pmx_saliency_map[k](j,i) = (float)rf.pdf(k);
			}
}



/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
MCMC_Particle_SurfaceColor::MCMC_Particle_SurfaceColor()
{
	// in L, u and v spaces
	for (int i=0; i<3; i++)
	{
		p_surface_model[i] = NULL;
		p_synthesized[i]   = NULL;
	}
}

MCMC_Particle_SurfaceColor::~MCMC_Particle_SurfaceColor()
{
	// in L, u, and v spaces
	for (int i=0; i<3; i++)
	{
		if (p_surface_model[i] != NULL)
			delete [](p_surface_model[i]);
		p_surface_model[i] = NULL;

		if (p_synthesized[i] != NULL)
			delete [](p_synthesized[i]);
		p_synthesized[i] = NULL;
	}
}

void MCMC_Particle_SurfaceColor::FitSurface(void)
{
	int			m;

	for (m=0; m<3; m++)
	{
		if (p_surface_model[m] != NULL)
			delete [](p_surface_model[m]);
		p_surface_model[m] = NULL;
		p_surface_model[m] = new MCMC_ControlPoints[num_particle];
		// free synthesized images
		if (p_synthesized[m] != NULL)
			delete [](p_synthesized[m]);
		p_synthesized[m] = NULL;
		// allocate space
		p_synthesized[m] = new Raster<float>[num_particle];
	}

	// original data 
	Raster<float> *pras[3];
	
	pras[0] = &(prespFilters->ras_L);
	pras[1] = &(prespFilters->ras_U);
	pras[2] = &(prespFilters->ras_V);

	for (m=0; m<3; m++)
		for (int k=0; k<num_particle; k++)
		{
			p_surface_model[m][k].Initialize(mx_map_temp,k,NULL,*(pras[m]));
			p_surface_model[m][k].Estimate();
		
			p_synthesized[m][k].SetDimension(prespFilters->rows(),prespFilters->cols());
			for (int j=0; j<prespFilters->rows(); j++)
				for (int i=0; i<prespFilters->cols(); i++)
					p_synthesized[m][k](j,i) = (float)p_surface_model[m][k].SurfaceValue(j,i);
		}
}


void MCMC_Particle_SurfaceColor::Set(MCMC_FilterResp &respFilters, mcmcMatrix<int> &clst, char *pfile_name)
{
	// initialization
	prespFilters = &respFilters;

	mx_map_temp = clst;
	num_particle = mx_map_temp.GetMax()+1;
	
	for (int m=0; m<3; m++)
	{
		if (p_surface_model[m] != NULL)
			delete []p_surface_model[m];
		p_surface_model[m] = new MCMC_ControlPoints[num_particle];
	}
	
	FitSurface();

	EM();

	CreateSalincyMap();
}

void MCMC_Particle_SurfaceColor::EM(void)
{
	int		i,j,k;
	double	dtemp;
	RF		rf;

	rf.Set(0,num_particle-1,1);

	for (int step=0; step<3; step++)
	{
		for (j=0; j<prespFilters->rows(); j++)
			for (i=0; i<prespFilters->cols(); i++)
			{
				for (k=0; k<num_particle; k++)
				{
					dtemp = DistanceToAParticle(prespFilters->ras_L(j,i),prespFilters->ras_U(j,i),prespFilters->ras_V(j,i),k,j,i);
					if (dtemp>1e-6)
						rf.SetOneBin(k,1/dtemp);
					else
						rf.SetOneBin(k,1e6);
				}
				k = (int)rf.SamplingNonSetCdf();
				mx_map_temp(j,i) = k;
			}
		// setp 2:
		FitSurface();
	}

	// any two clusters are overlapping?
	mcmcVector<int> vc_sign(num_particle);
	vc_sign.InitValue(1);
	for (i=0; i<num_particle; i++)
	{
		if (vc_sign(i)==1)
		{
			for (j=i+1; j<num_particle; j++)
			{
				if (vc_sign(j)==1)
				{
					if (DistanceBtwnParticles(i,j)<3.0)
						vc_sign(j) = 0;		// tow clusters are too close
				}
			}
		}
	}
	mcmcMatrix<int> vc_index(num_particle,1);
	double		dist_temp;
	int			ipar_temp,m;

	m = 0;
	for (i=0; i<num_particle; i++)
		if (vc_sign(i)==1)
			vc_index(i) = m++;
	
	// relabel
	for (j=0; j<prespFilters->rows(); j++)
		for (i=0; i<prespFilters->cols(); i++)
		{
			dist_temp = 1e90;
			ipar_temp = -1;

			for (k=0; k<num_particle; k++)
			{
				if (vc_sign(k)==1)
				{
					dtemp = DistanceToAParticle(prespFilters->ras_L(j,i),prespFilters->ras_U(j,i),prespFilters->ras_V(j,i),k,j,i);
					if (dtemp<dist_temp)
					{
						ipar_temp = k;
						dist_temp = dtemp;
					}
				}
			}
			assert(ipar_temp>=0);
			mx_map_temp(j,i) = vc_index(ipar_temp);
		}

	num_particle = m;

	FitSurface();


	Raster_Color<uchar> aa;
	char		  name[200];
	aa.SetDimension(prespFilters->rows(), prespFilters->cols());
	for (k=0; k<num_particle; k++)
	{
		for (j=0; j<prespFilters->rows(); j++)
			for (i=0; i<prespFilters->cols(); i++)
				aa.SetPixelInLUV(j,i, (int)(p_synthesized[0][k](j,i)), 
						         (int)(p_synthesized[1][k](j,i)),(int)(p_synthesized[2][k](j,i)));
		sprintf(name,"a_%d.ppm",k);
		WriteToPPM(aa,name);
	}
}

double MCMC_Particle_SurfaceColor::DistanceToAParticle(const double pixel_l, const double pixel_u, const double pixel_v,
													   const int ipar, const int j, const int i)
{
	double dl,du,dv;
	dl = p_synthesized[0][ipar](j,i);
	du = p_synthesized[1][ipar](j,i);
	dv = p_synthesized[2][ipar](j,i);

	return sqrt(pow(dl-pixel_l,2.0)+pow(du-pixel_u,2.0)+pow(dv-pixel_v,2.0));
}

double MCMC_Particle_SurfaceColor::DistanceBtwnParticles(const int ipar1, const int ipar2)
{
	double	dis1,dis2;
	int	num1,num2;
	double	dl1,du1,dv1,dl2,du2,dv2,dtemp;

	//initialization
	dis1 = 0.0;
	dis2 = 0.0;
	num1 = 0;
	num2 = 0;

	for (int j=0; j<prespFilters->rows(); j++)
		for (int i=0; i<prespFilters->cols(); i++)
		{
			dl1 = p_synthesized[0][ipar1](j,i);
			du1 = p_synthesized[1][ipar1](j,i);
			dv1 = p_synthesized[2][ipar1](j,i);

			dl2 = p_synthesized[0][ipar2](j,i);
			du2 = p_synthesized[1][ipar2](j,i);
			dv2 = p_synthesized[2][ipar2](j,i);
			
			dtemp = sqrt(pow(dl1-dl2,2.0)+pow(dl1-dl2,2.0)+pow(dl1-dl2,2.0));
			if (mx_map_temp(j,i)==ipar1)
			{
				dis1 += dtemp;
				num1++;
			}
			else if (mx_map_temp(j,i)==ipar2)
			{
				dis2 += dtemp;
				num2++;
			}
		}
	if (num1 <= 0)
		dis1 = 1e5;
	else
		dis1 = (dis1/(double)num1);
	if (num2 <= 0)
		dis2 = 1e5;
	else
		dis2 = (dis2/(double)num2);
	return min(dis1,dis2);
}


void MCMC_Particle_SurfaceColor::CreateSalincyMap(void)
{
	RF		rf;
	int		i,j,k;
	double	dtemp;

	if (pmx_saliency_map != NULL)
		delete []pmx_saliency_map;
	pmx_saliency_map =  new mcmcMatrix<float>[num_particle];
	for (k=0; k<num_particle; k++)
		pmx_saliency_map[k].SetDimension(prespFilters->rows(),prespFilters->cols());

	rf.Set(0,num_particle-1,1);
	
	for (j=0; j<prespFilters->rows(); j++)
		for (i=0; i<prespFilters->cols(); i++)
			{
				for (k=0; k<num_particle; k++)
				{
					dtemp = DistanceToAParticle(prespFilters->ras_L(j,i),prespFilters->ras_U(j,i),prespFilters->ras_V(j,i),k,j,i);
					if (dtemp>1e-6)
						rf.SetOneBin(k,1/dtemp);
					else
						rf.SetOneBin(k,1e6);
				}
				rf.normalize();
				for (k=0; k<num_particle; k++)
					pmx_saliency_map[k](j,i) = (float)rf.pdf(k);
			}
}
